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We propose a comprehensive dynamical model for cooperative motion of self-propelled particles, 
e.g., nocking, by combining well-known elements such as velocity-alignment interactions, spatial 
interactions, and angular noise into a unified Lagrangian treatment. Noise enters into our model 
in an especially realistic way: it incorporates correlations, is highly nonlinear, and it leads to a 
unique collective behavior. Our results show distinct stability regions and an apparent change in 
the nature of one class of noise-induced phase transitions, with respect to the mean velocity of the 
group, as the range of the velocity-alignment interaction increases. This phase-transition change 
comes accompanied with drastic modifications of the microscopic dynamics, from nonintermittent to 
intermittent. Our results facilitate the understanding of the origin of the phase transitions present 
in other treatments. 
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I. INTRODUCTION 



The study of the collective motion that emerges from 
short-range interactions in systems of self-propelled par- 
ticles (SPPs) is of great interest nowadays due to po- 
tential applications in physics, engineering, and biology 
USUI- In nature, the formation of bird flocks or animal 
herds implies the occurrence of a condensation in veloc- 
ity space and a condensation in position space: the con- 
stituents of the group move with velocities similar to one 
another and additionally, tend to form a spatially con- 
tiguous collection. Such behaviors have been classified as 
velocity matching or alignment, flock centering, and colli- 
sion avoidance or separation [J] . The additional element 
of noise is always present in realistic systems, nonethe- 
less, the way it is introduced may have non-trivial effects 
in the modeling of flocking phenomena [H, Models 
in the literature have concentrated in all or subsets of 
these ingredients @, & 1 EE El El 0, El EE El, E3 • 
However, they have not studied the effects of nonlinear 
correlated noise, which in the model presented in this 
article appears in a natural way. 

The physics of dynamic phase transitions is ideally 
poised to address cooperative phenomena, hence, the 
considerable activity in this field in recent times [l8l |. 
One may argue that flocking of real biological systems 
in nature occurs in continuous time, that alignment is 
highly important to the process, and that complex non- 
linear noise, even correlated, is an essential ingredient 
of realistic systems. This argument is the motivation 
for our study. Succinctly stated, in this article we study 
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the effects of alignment interactions, centering and corre- 
lated angular noise, in a Lagrangian (trajectory-based) 
description of flocking, meaning that dynamical equa- 
tions of motion are employed to compute the trajectories 
of the individual particles. By doing so, we provide a 
rich description of the different dynamics of the system, 
including the origin of intermittent behavior displayed by 
the mean velocity of the group. 

This article is organized as follows. In Sec. [Til we intro- 
duce the model of our study. We explain how noise enters 
naturally non-linear and correlated, as we describe how 
the alignment interaction is implemented. In Sec. IIIII wc 
define the quantities we use to characterize the system, 
and report the results as we discuss the origin of the dif- 
ferent dynamical regimes observed. Finally, in Sec. IIVI 
we compare our results in qualitative terms with some 
other models of interest and report our conclusions. 



II. MODEL 

We study N interacting particles, each having position 
vector x.; and velocity vector v,; = dx.i/dt obeying 
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where m and vq are the mass and the preferred speed of 
each particle respectively, K is the coupling constant for 
the real-space condensation interaction, and 7 is the am- 
plitude of the self-propulsion term with propulsive direc- 
tion ipi and noise 4>i . This term derives from the assump- 
tion of constant (resultant) propulsion force of Czirok et 
al. [ill ], later generalized by Levine et al. 13] to incorpo- 
rate in a straightforward manner the alignment interac- 
tion by directing this force in the motion's mean direc- 
tion ipi of the neighbors of a given particle i. We modify 
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this term further in order to introduce correlated angular 
noise. This is achieved by simply adding the stochastic 
variable <pi(t) to ipi, as shown in Eg. ([1]), much in the 
spirit of (but not equivalent to) Ref. 0- 

In order to determine the mean direction ipi of the 
neighbors of a given particle i, the alignment interaction 
is implemented by allowing particle i to interact only 
with its /i nearest neighbors at time t, thus, defining the 
set Q(t) = {v n (7j),v. i2 (/j), . . . ,v i(j (t)} (v t Q) with the 
velocity vectors of the neighbors. We identify /i as the 
connectivity of the particles, and consider it to be equal 
for all particles. In this way, given that the individual 
velocities have a magnitude Vi and direction ipi, the 
case fj, = corresponds to the absence of alignment in- 
teractions, and ipi = ipi, otherwise 



ipi = arctan 



Vi sin 
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On the other hand, noise introduced through the 
stochastic variable ipi , as in Eq. |T]) , is intrinsic to the de- 
cision mechanism of each one of the self-propelled agents 
in that a given agent, receives a clear signal from its 
neighbors, however, it may "decide" to move in a differ- 
ent direction Q. For example, this kind of noise may be 
thought as arising from the limitations of a bird to fol- 
low a specific and precise flight-path (due to fatigue, for 
example) or some other external factors (wind currents, 
for instance). In our model, noise appears naturally in 
a highly non-linear form requiring correlations (that de- 
cay in time) to be considered for it to have a non-trivial 
effect [19] j_and differs fundamentally from previous stud- 
ies @,TE S, EI El El EE El H|- In this work, the 
correlated random sequences <pi(t), yet independent be- 
tween different agents, are obtained by correspondingly 
coupling the stochastic equation 



K =-X<f>i + Ag t , 



(3) 



to each one of the equations of motion of the group in- 
dividuals. Here, gi(t) is a Gaussian white noise with 
mean value (gi(t)) = and autocorrelation function 
(gi(t)gj(t')} =2D6ij5(t-t'). Thus, the driven noise 
(p'i{t) is then exponentially correlated noise (derived 
from the Orstcin-Uhlcnbcck processes) with properties 
(#(*)) = and {{ti(tW j (t))} = D\5 ij eiq ? (-\\t-t>\). 
The curly brackets { • • • } denote averaging over the dis- 
tribution of initial (p' Q values, which are equivalent for 
every particle, and taken from the distribution 
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to compare our results to some other models of interest, 
we would like for our correlated angular noise, (pi{t), to 
have an uniform distribution 



k if \<P\<V 
otherwise 
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in the interval [—77, 77], where rj (running from to 7r) de- 
fines the maximum possible value for \(pi\. By matching 
the mean squared variance of the distribution of Eq. 
(j4|) with that one of the desired distribution P v of Eq. |5]), 
i.e., DX = rf /3, with a further transformation that pre- 
serves the area between them, the correlated sequences 
of angular noise can be obtained from 



(6) 



where erf(x) stands for the Gauss error function. For a 
given (dimensionless) noise intensity 77, and correlation 
time t c , the properties of our angular noise are 
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which correspond to its mean value and autocorrelation 
function respectively. 

Throughout the rest of this paper, we will measure 
speed in units of vq, time in units of to = mvo/j and 
distance in units of xq = ttivq/j. In this way, only three 
parameters remain independent: the dimensionless cou- 
pling parameter k = itiRvq/j 2 , the connectivity /1, and 
the intensity of the angular noise 77. All numerical re- 
sults presented here were obtained by integrating the set 
of Eqs. (H|) and ^ with a modified version of the velocity- 
Verlet algorithm (22J. 



III. DISCUSSION AND RESULTS 

To characterize the system and its stationary states, 
we monitored quantities such as the instantaneous po- 
sition and velocity of the center of mass, defined as 

X C m(*) = iEL x «( ( ) and V C m(*) = ^Ejli v i(*)> 
respectively. To quantify the translational and rotational 
motions of the flock, as order parameters, we measured 
the normalized mean velocity and angular momentum of 
the group, 
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Clearly, the correlation time for the colored noise is 
t c = A -1 . Notice that the noise sequence has a 

Gaussian distribution [21] . However, in order to be able 



where l>i(t) — Xj(t) x Vj(t) corresponds to the angu- 
lar momenta of the individual particles with respect 
to the center of mass, Xj(i) = x,(t) — XcAf(i), and 
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FIG. 1: (Color online) Snapshots of the configuration in posi- 
tion space for the model given in Eqs. ^ and ((2)1 , with k = 1 
and t c = 10, for systems of N = 100 particles. The small 
arrows represent the velocity vectors Vi for each one of them 
while the big (red) arrows, across figures, correspond to the 
way of increasing noise intensity. In (a) and (b), and (c) and 
(f), the translational and oscillatory states are, respectively, 
shown in regions I and IV for different connectivity values (see 
text for more details). When the system shows intermittent 
dynamics, (d) a laminar regime can be interrupted by (e) a 
turbulence burst. In (d) and (e), the value of n considered 
was taken inside region II (see text). 
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FIG. 2: (Color online) Stationary order parameter ^, given 
in Eq. (|8a[l . as a function of the noise intensity r\. Equations 
|T|) and ([2]) were numerically solved for systems of N — 200 
particles with re = 1 and t c = 10. Initial conditions were taken 
in the TranS (a) and in the OscS (b). The changes of stability 
and the different kinds of phase transitions are apparent as 
the connectivity fi/N changes. See text for more details 



v;(t) = Vi(i) — VcAi(t). We also measured quantities 
regarding the spatial distribution of the particles such as 
the instantaneous mean square dispersions in the direc- 
tions parallel (£||) and orthogonal to the direction 

0fV C M H, 

= N{ ^ Mit)l2 hm^CM(t)r, (9a) 

1 N 

E±(t) = N\V CM (tW ^ [Mt)XVcM{t)]2 - (9b) 

Depending on initial conditions, the values of the nor- 
malized connectivity fi/N and the noise intensity ?y, the 
system may attain a translational state (TranS) where all 
of the particles move more or less in the same direction, 
shown in Figs. [Ha) andQJc), with >F > 0. The curved 
shape of the group is principally due to correlations in the 
noise in combination with the centering interaction, for 
it can be observed even in the absence of any alignment 
interactions. The system may also attain an oscillatory 
state (OscS) shown in Figs.QJb) andQJf), where particles 
circle around a center of mass that diffuses with $«0. 

Particular care was taken to obtain initial conditions 
for the system either in the TranS or the OscS by setting 
7i = at the beginning. Then, in one case, the system 
was allowed to relax to a stationary TranS under a global 
alignment interaction condition fi/N = N — 1. In the 
other case, particles were homogeneously placed on a ring 
with tangent velocity vectors, and the system was allowed 
to relax to a stationary OscS with n/N = 0. At this 



point, time is set t — in all our numerical simulations, 
and the parameters of the system are correspondingly 
fixed for each one of the cases studied. 



A. Mean velocity of the group 

In the absence of alignment interactions (fJ,/N = 0), 
the normalized mean velocity of the group, 'F, undergoes 
a noise-induced discontinuous phase transition from the 
TranS to the OscS as shown in Fig. [5] where plots for 'F 
as a function 77 are presented. Both the TranS and the 
OscS are accessible and stable for subcritical values of the 
noise, as can be appreciated in Figs.EJa) and[2{b), where 
initial conditions were taken in the TranS and the OscS, 
respectively; this kind of behavior resembles the model 
of Erdmann et al. 14]. The critical point r] c w 0.548 
of the phase transition for the case shown in the figure, 
shifts to lower values as one increases t c due to stronger 
effects from the noise. In contrast, increasing the cou- 
pling parameter with the center of mass, re, enhances the 
stability of the TranS and r/ c shifts to higher values. On 
the other hand, when a global alignment interaction is 
considered, i.e., when fi = N — 1, the TranS solution be- 
comes the only available solution of the system as the 
phase transition disappears. For this case, it can be an- 
alytically estimated that $ w ^p; we analyze this and 
some other limit cases in Appendix [X] Nevertheless, for 
local alignment interactions with fi/N <C 1, at least four 
regions with different stability can be identified, depicted 
with roman numerals in Fig. [21 
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FIG. 3: (Color online) Times series of the order parameter 
= |VcA/(i)|/«o for different values fi/N and rj. Typical 
cases inside region III for different connectivities are shown in 
(d), where the curves are vertically displaced for better clarity. 
Note that the time intervals between turbulence bursts, and 
between free transitions to the translational and oscillatory 
states inside reigon III, are in general larger than the time 
scale of the noise-induced fluctuations for fi/N < 0.2. See the 
text for more details. 



Region I (gray checked portions in Fig. [2] for the 
case /i/N = 0.02) shows bistability similar to the case 
fi/N = 0, where both states are accessible depending 
on initial conditions. For larger values of ij, we found 
a window of noise values (region II) where the TranS 
is the only available stable solution, whose size increases 
with the normalized connectivity fi/N. Indeed, the noise- 
induced fluctuations make the OscS unstable inside this 
region even if initial condition are taken in this state, 
consequently driving the system to the TranS [see, for 
example, the event marked with the downward arrow for 
the curve with -q — 1.45 in Fig. |3fa)]. Thus, the system 



TABLE I: Approximate values of the critical points that sepa- 
rate the different stability regions of the cases with finite fi/N 
analyzed in Fig. [2ja) . 



fi/N r/ot rj tb 
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FIG. 4: (Color online) Stationary probability distribution 
functions of the order parameter <& for selected cases of Fig. [2] 
with finite fi/N. In (a)-(c), the filled curves correspond to val- 
ues of r\ inside region III, where the system shows bistability 
and free transitions between the TranS and the OscS through 
the bimodal distribution of P(^\ This behavior is typical 
of first-order phase transitions but only close to the critical 
point. In (d), the filled curve denotes the unimodal character 
of P(\&) in the transition from the TranS to the OscS, typical 
of second-order phase transitions. 



exhibits a discontinuous phase transition in , F from re- 
gion I to region II with critical point r] t; the values of 
the critical points for the different cases are presented in 
Table D 

Inside region II, the solution for the TranS starts to 
shows disorder spikes as the system develops intermittent 
behavior [some of these events are marked with upward 
arrows for the curve rj = 1.611 in Fig. [3[a)]. Increasing 
the connectivity and/or the noise results in an increase in 
the amplitude of the spikes and, in some cases, the system 
is able to even segregate into clusters. See, for example, 
how for a fixed and the same noise intensity, the laminar 
regime of the TranS shown in Fig. [ljd) is interrupted by 
a segregation event shown in Fig. [lje) . These events are 
reminiscent of "turbulence bursts" where the structure 
of the group is directionally inhomogeneous; however, if 
present, the clusters themselves may show local order [see 
Fig. [He)]. With a further increase in the noise intensity, 



and for values up to fi/N w 0.2, a region with mixed dy- 
namics develops in a discontinuous manner (region III), 
depicted by the shadowed portions in Fig. [2] for the cases 
n/N = 0.02,0.1,0.2. As can be appreciated, the size 
of this region increases with the normalized connectivity 
fx/N. 

Inside region III, the stronger turbulence bursts may 
cause the system to suffer free transitions between the 
TranS to the OscS, in the sense that the system may 
acquire any of the two states regardless of initial condi- 
tions as shown in Fig. [3Jb) for the curve with r\ = 2.094. 
We must mention that turbulence bursts in systems that 
present intermittent behavior are known to occur at ir- 
regular time intervals showing a power-law distribution 
[20l [23l [24j and, in our model, not all of the turbulence 
bursts induce a transition between the TranS and the 
OscS. What we observed from our numerical results is 
that, as the connectivity increases, the frequency of these 
transitions also increases, while the plateaus where the 
system spends some time in either of the two states be- 
come smaller in the mean [see Fig. Eld)]. Moreover, the 
system spends more time in the TranS at the beginning 
of region III, while, with the increasing noise, it starts to 
spend more and more time in the OscS as one approaches 
the end of this region. This can be better appreciated 
looking at the probability distribution function (PDF) 
of * shown in Figs. Ha) and g^b) for p/N = 0.02,0.1. 
There, the bimodal filled curves correspond to values of 
t] inside region III, and corroborate the coexistence of 
states and the bistability. The rightmost peaks of the 
distributions correspond to the TranS, and their ampli- 
tude diminishes with the increasing noise. The opposite 
happens to the leftmost peaks that correspond to the 
OscS, whose amplitude increases with the noise intensity. 
Nonetheless, this "transference of stability" between the 
translational and oscillatory states occurs in a continu- 
ous manner as shown in Figs. [Ha) and[^b); the strong 
fluctuations for the cases with fi/N — 0.02,0.1 are due 
to lack of better statistics for the free transitions, requir- 
ing longer integration times than those considered in this 
article. 

It is not our intention to give a quantitative detailed 
analysis of the distribution of the time intervals between 
the free transitions inside region III here due to the cum- 
bersome computational requirements for these calcula- 
tions. We will leave this as an open question for future 
work. 

The borders of region III are delimited by the criti- 
cal points rjtb, for the phase transition that occurs from 
regions II to III, and rjbo for the one that occurs from 
regions III to IV. The values of the critical points for the 
different cases presented in Fig.[5]are shown in TableU In 
particular, the end of region III corresponds to the case 
when the OscS becomes the only stable solution of the 
system in region IV, in a contrasting effect from the cor- 
related angular noise in comparison to region II. These 
transitions from regions II to III, and III to IV, may be 
considered discontinuous as the coexistence of states in- 
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FIG. 5: (Color online) Stationary angular momentum A, 
given in Eq. ()8b|) . as a function of n with N = 200, K = 1, 
and t c = 10, for initial conditions taken in the TranS (a) and 
in the OscS (b). For the case fx/N = 0, the vertical dotted 
line corresponds to the critical point r\ c of the phase transi- 
tion in $ from the TranS to the OscS when initial conditions 
are taken in the TranS. Beyond that line the only solution 
of the system is the OscS regardless of initial conditions. On 
the other hand, beyond the vertical dash-dot-dashed line, the 
correlated angular noise is strong enough to randomly change 
the direction of rotation of the particles in the OscS. In con- 
trast, in between the two vertical lines and after the system 
reaches the stationary OscS, particles rotating in either di- 
rection, clockwise or counterclockwise, will remain rotating 
in such a way regardless of the fluctuations induced by the 
noise. This is more evident in (a) where A presents some 
bumps product of an asymmetric number of particles rotating 
in both directions, due to the transition the system undergoes 
from the TranS to the OscS. See text for more details on the 
cases with finite connectivities. 



side region III is reflected in the bimodal distribution of 
P(^) as explained before; see, for example, Figs. Ufa) - 
2{c). Nonetheless, this contrasts to what typical first- 
order phase transitions show, where the bistability and 
the coexistence of states can only be observed close to the 
critical point and not for a continuous window of values 
of the noise. 

For connectivities fx/N > 0.3, a new kind of dynamics 
develops as the frequency of the occurrence of turbulence 
bursts and free transitions between the translational and 
oscillatory states become indistinguishable from the time 
scale of the fluctuations induced by noise. For this case, 
^> smoothes out completely and, as region III disappears, 
the transition from the TranS in region II to the OscS 
in region IV shows an apparent continuous nature with 
critical point r) to (see Fig. [2] and Table UJ). The nature of 
this phase transition can be inferred from the one-peaked 
form of P(^f) at criticality shown in Fig.UJd), and from 
the fact that ^(t) does not show clear plateaus that im- 
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FIG. 6: (Color online) Probability distribution functions of 
the magnitude of the angular momentum of the individual 
particles Li for systems with N = 200. In (a), for fi/N = 0, 
the two peaks away from the center for the curve with r\ — 
0.564 indicate a limit-cycle oscillatory state. That is also the 
case in (b), for the different values of n/N with a one-peaked 
P\(L) away from the center. In (c) and (d), the filled curves 
correspond to values of r\ inside region III. In (c) the tails in 
some of the distributions indicate the development of limit- 
cycle- like oscillatory states. In contrast, in (d), the tails of 
the distributions decay faster from the center indicating the 
development of qualitatively disorder oscillatory states. 



ply the coexistence of states but only strong fluctuations 
[see the curve for 77 = 2.175 in Fig. |3Jc)] as the system 
goes from the TranS to the OscS. Looking at the angu- 
lar momentum of the group can provide a further insight 
(see next subsection). 

Numerical results for different values of k (not shown) 
indicate that the system displays the same qualitative 
behavior, and the same stability regions, just as in the 
case considered here with k = 1. Moreover, the critical 
points that separate these regions behave similarly to the 
critical point ij c for the case fi/N = 0, i.e., they shift 
to higher values as k increases while they shift to lower 
values as K decreases. Nonetheless, it was also observed 
that lager values of the connectivity fx/N are required 
for the phase transition from regions II to IV to become 
apparently continuous as k increases. 



B. Angular momentum 

Figures [SJa) and M,b) show plots of the angular mo- 
mentum A as a function of 77, with initial conditions taken 
in the TranS and the OscS, respectively. For fi/N = 0, A 
fluctuates around zero for stationary translational states. 
On the other hand, for stationary oscillatory states, par- 
ticles oscillate around the center of mass, on a limit-cycle, 
much in the same way as oscillatory states present in 
some other models [H GJ, Q GJ, (TJ • This can be ap- 
preciated in Fig. [nja) from the two peaks displayed by 
the stationary PDF, P\(L), of the magnitude of the in- 
dividual angular momenta Li = |Lj| of the particles for 
77 = 0.564. 

Having in mind that the OscS is always a solution for 
the system when fi/N = 0, becoming the only solution 
for 77 > rj c , the value of the order parameter A will de- 
pend on initial conditions for noise values (77 < 0.867) 
to the left of the dash-dot-dashed vertical line in Figs. 
EJa) and [5£b), whenever a stationary oscillatory state 
is reached. This can be understood from the fact that, 
in this region, and in the absence of alignment interac- 
tions, the correlated angular noise is not strong enough 
to randomly change the particles direction of rotation. 
In particular, in Fig.[5fb), initial conditions for the OscS 
were taken on the limit-cycle with half of the particles 
rotating in the clockwise direction, and half in opposite 
one, resulting in a A — 0. Nonetheless, when initial con- 
ditions are taken in the TranS, but for noise values 77 > r\ c 
[to the right of the dotted vertical line in Fig. Ufa)], A 
presents some bumps. Under these conditions, for val- 
ues of the noise (rj c < rj < 0.867) between the dotted 
and dash-dot-dashed vertical lines in Fig. [S£a), the sys- 
tem undergoes a transition to the OscS after spending a 
transient time in the TranS (also reported in Ref. 
Thus, once the system reaches a stationary OscS, the 
number of particles rotating in each direction (clockwise 
and counterclockwise) may be asymmetric deriving in a 
fixed A > 0. In contrast, for noise values (77 > 0.867) to 
the right of the dash-dot-dashed vertical line, where the 
correlated angular noise is strong enough to randomly 
change the particles direction of rotation in the OscS, 
A«0 regardless of initial conditions. However, for finite 
connectivities (n/N > 0), A shows the different stability 
regions identified before. 

In region I, as the alignment interactions induce the 
majority of the particles to randomly select a direction 
of rotation, Awl for the OscS. This state has a limit- 
cycle quality as P\{L) shows in Fig. [S^b). In region II 
only the TranS is accessible, thus, the value of A is closer 
to zero. 

Inside region III, and for connectivities up to fJ,/N rs 
0.2, A shows a discontinuous peak as the majority of the 
particles tend again to rotate in the same direction when- 
ever the system attains an OscS. This is corroborated by 
the tails of P X (L) in Fig. E^c) for fi/N = 0.02; however, 
the direction of rotation is selected randomly considering 
that every OscS follows a free transition from the TranS. 
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FIG. 7: (Color online) Above, (a)-(c), mean-square disper- 
sions, given in Eqs. (|9a|l and (|9b[) . and the mean range of the 
alignment interaction, (ro) (see text for definition), as a func- 
tion of rj, for different connectivities, N = 200 and t c = 10. 
The curves with dashed lines and solid symbols correspond to 
numerical results with initial condition in the OscS, while the 
curves with solid lines and clear symbols correspond to numer- 
ical results with initial condition in the TranS. The different 
stability regions are identified with roman numerals for the 
different cases. Below, behavior of the maximum transversal 
dispersion, Ef AX , as a function of (d) t c , and (e) ji/N. 



This type of OscS resembles a limit-cycle one even though 
the tails in P\(L) are wider than in the previous cases 
due to particles that transit close to the center of mass 
[see Fig.QJb)]. The presence of a discontinuous peak in A 
corresponds to cases where P(\&) shows a bimodal distri- 
bution inside region III. Nonetheless, with the increasing 
connectivity, the amplitude of this peak decreases (see 
Fig. [SJ along with the "limit-cycle" quality of the OscS. 

Of particular interest is the case with jx/N = 0.2 since, 
inside region III, it presents the discontinuous peak in A 
along with a bimodal distribution for P(^S); however, the 
distribution for P\(L) looks rather symmetric through- 
out regions II to IV as shown in Fig. [Hd). We will use 
this fact to define a crossover regime that separates the 
cases when the transition from region II to IV presents 
the bistable region III, from those when it seems contin- 
uous and region III is absent (see next subsection). 
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FIG. 8: (Color online) Snapshots of the configuration in posi- 
tion space for the model given in Eqs. (JTJ and ^ with k = 1, 
t c — 10, and N — 200. The small arrows represent the veloc- 
ity vectors Vj for each one of the particles. These states were 
obtained from simulations with initial conditions in the OscS. 
In (b) the large arrow points to the trajectory of the center 
of mass (solid line) for a finite period of time. 

In region IV, the values of A are closer to zero as the 
OscS becomes qualitatively disordered due to stronger 
noise, while P\(L) does not show any tails but a symmet- 
rical distribution around zero, in particular for ji/N > 
0.1. For fi/N > 0.3, as region III disappears and the 
transition between the TranS and the OscS becomes ap- 
parently continuous, the discontinuous peak in A also dis- 
appears, giving way to a smooth dependence of A on rj 
instead. In this case, the OscS does not show any limit- 
cycle quality, looking rather disordered, while P\(L) is 
always symmetrical throughout the transition, similar to 
the case for ji/N = 0.2. Thus, the apparently continuous 
transition from the TranS to the OscS, for connectivities 
n/N > 0.3, may be considered an order-disorder. 



C. Size of the system, dispersions, and critical 
points. 

In addition, the values of t c used in this work corre- 
spond to a regime where the size of the flock in position 
space, for any given ry, does not depend on t c anymore. 
In consequence, the dependence of the critical points that 
separate the different stability regions on t c becomes also 
constant for a fixed values of k and ji/N. This regime 
is defined by the condition t c ^> 1, and can be explained 
by looking at the behavior of dispersions ©, shown in 
Figs. E[a) - [7£c) as a function of ?/. In particular, by 
looking at the maximum of the transversal dispersion, 
S" AX , and its dependence on t c ; as shown in Fig. [3d), 
it becomes constant for t c ^S> 1. In the same limit, £" AX 
does not change much with ji/N but rather lightly with 
N [see Fig.JTJe)]. This means that the size of the system 
depends mainly on the noise intensity rj and the corre- 
lation time t c , in the end becoming independent on t c 
for t c ^ 1. Thus, the critical points that separate the 
different stability regions will mainly depend on the ra- 
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FIG. 9: Plot fj,/N vs. TV where the solid line with grey circles 
corresponds to the crossover regime that separates the cases 
where the system exhibits a bistable region (region III) with 
free transitions between the TranS and the Oscs (depicted 
with clear circles), from those where the transition from the 
TranS to the OscS is apparently continuous (depicted with 
solid circles). In all of the cases t c = 10 and k = 1. 



tio /i/N since the mean range of the alignment interac- 
tion, (r D )(t) = iEZii^E^ecmlMt) ~ *At)\], still 
depends on rj and fi/N itself, as can be appreciated in 
Figs.EJa) -E{c), and not on t c (in general, (r a ) is a func- 
tion of all the parameters of the system). This is one 
of the main effects of the correlated angular noise, and 
allows an increase in the normalized connectivity fi/N 
to translate into an effective increase in the range of the 
alignment interaction, in position space, with respect to 
the overall size of the group for any given value of r\ as 
shown in Figs. Eta) -[2c). 

We must remark that for cases with larger connectiv- 
ities fJ>/N > 0.2, when initial conditions are taken in the 
OscS and for weak noises (rj <C 1), the resulting station- 
ary oscillatory states can present the formation of clusters 
as shown in Fig. [8] This comes as a result of the short- 
range character of the alignment interaction. Nonethe- 
less, structures like these have been obtained by provid- 
ing each particle with a short-range repulsive potential 
|15j . On the other hand, the clusters, when present, 
seem to be very stable since the exchange of particles 
is not likely as they oscillate around the center of mass 
in a limit-cycle fashion. The segregation of the system 
in this way is the main cause that the dispersions par- 
allel (E||) and orthogonal given in Eqs. ©, are 
different in Figs. Et D ) andEtc) for weak noise intensities. 
What happens is that, since the particles are not homo- 
geneously distributed in a circular form around the center 
of mass, the center of mass itself does not lie in the center 
of the circular formation but closer to the largest cluster. 
In consequence, the center of mass is able to perform a 
translational motion rather than a diffusive one. This is 
shown in Fig. [5]Jb) , where the large arrow points to the 
trajectory of the center of mass for a finite time. The ex- 
change of particles among the clusters is more likely when 
noise increases, in the end leading to the destruction of 
these structures. 

Finally, Figure [H] shows a plot of fi/N vs. N, and is in- 



tended as a scaling analysis of the transition from regions 
II to IV. One can observe that as the number of particles 
N in the system increases, the crossover regime (solid 
line with grey circles) that separates the cases when the 
bistable region III is absent, and the transition from the 
TranS to the OscS seems to be continuous with a one- 
peaked P(\&) at criticality (solid black circles above the 
line), from those when a region III is present and P{9) is 
bimodal (clear circles below the line), moves from lower 
values of the connectivity fi/N to larger ones, finally sta- 
bilizing for N > 200. As far as our numerical results 
indicate, increasing N any further docs not change the 
position of the different stability regions for finite values 
of n/N and t c > 1. 



IV. CONCLUSIONS 

We have studied the mixed effects of correlated angu- 
lar noise, centering and alignment interactions in a sim- 
ple two-dimensional (2D) Lagrangian model of N self- 
propelled particles. As a result, we have observed rich 
dynamics in the collective behavior along with different 
stability regions as the amplitude of the noise changes. 
Of particular interest is the development of intermittent 
behavior as the presence of alignment interactions induce 
turbulence bursts in the TranS where the group is able 
to even segregate into clusters. During these events, the 
group may present local order but not global, and for 
certain combinations of the parameters, they allow the 
system to freely transit from the TranS to the OscS and 
vice versa. This result is distinctive from other models 
where the collective dynamics is condemned to live in 
monotonous stable states. Indeed, a closer inspection of 
P(^>) shows the presence of a tail that deviates from a 
Gaussian distribution and accounts for these turbulence 
bursts inside region II, while this behavior is not observed 
in the case fi/N = [see Figs. HUT a) - [TUTc)]. Moreover, 
increasing the range of the alignment interaction (by in- 
creasing n/N) enhances this effect that finally leads to 
a fundamental change in the dynamics of the system as 
the bistable region III disappears, and the transition in \ff 
from the TranS of region II to the OscS of region IV ap- 
parently changes its nature to continuous [see Figs.^a) 
and [S]. This comes about from the fact that, for large 
enough connectivities (/x/iV > 0.2) and subcritical values 
of r\ closer to the end of region II, the system becomes 
fully intermittent as P(^) shows a non-Gaussian asym- 
metrical form [see Fig. [TUTd)]. consequently leading to 
the existence of region III and its further disappearance 
when the connectivity increases. This strongly suggests 
what seems to be a clear connection between the inter- 
mittent behavior in the collective motion of the group 
and the kind of phase transitions undergone by its mean 
velocity. 

The kind of dynamics the system develops inside re- 
gion II, resembles that in the model of Vicsek et al. 0], 
where the cluster formation has proven to play an es- 
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FIG. 10: (Color online) [(a)-(d)] Semi-log plots of the PDF of 
* with N = 200, t c = 10, and k = 1, for different values fi/N 
and r\. In (a), the selected value of r] lies in the subcritical 
region of the phase transition from the TranS to the OscS. In 
(b)-(d), the selected values of r\ lie inside region II. The dotted 
lines correspond to Gaussian fits of the whole distribution in 
(a), and of the crest of the distributions in (b)-(d). 



sential role in the development of intermittent dynam- 
ics [13, S3]. This comes accompanied by a continuous 
order-disorder phase transition in the mean velocity of 
the group [IE SE EE S3, El ■ The intermittent behavior 
shows up as, in the stationary ordered phase, a laminar 
regime (where all the particles move in a rather ordered 
fashion) is interrupted by chaotic bursts of turbulence 
where the system segregates in clusters that may present 
local order, but not global, in a very similar way as our 
model does. Moreover, the authors in Ref. [13] reported 
the cluster formation in a Lagrangian model that con- 
siders alignment forces, centering, and separation, but 
with white additive noise. In fact, intermittent behavior 
can be observed in other systems far from equilibrium 
including biological ones (see references in [23jl. 

Finally, as our numerical simulations show, the dif- 
ferent dynamics displayed by our model seem to be ro- 



bust even if each particle is provided with a hard-core 
repulsive potential. This case is analyzed in detail in 
[29j. On the other hand, the nonlinear way angular 
noise enters in our model results in the development of 
an effective individual mean speed that depends on 77, 
and in contrast to some other models where it does not 
[13, El, EE E3, E3| ■ It can be shown that Vi fluctuates 
around the value — sin 77, since (u,) depends only on the 
form of the distribution of the noise. We believe our re- 
sults may be relevant in the theory of flocking and, in 
general, in the theory of phase transitions in systems out 
of equilibrium. 
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APPENDIX A: LIMIT CASES 

In this appendix we present the analysis of some limit 
cases of our model, and of other well-known models, as we 
compare them in the light of some analytical estimations. 

First, let us consider the global alignment interaction 
limit of model {l}, where the individual propulsive direc- 
tion ipi of each particle, given in Eq. ([2]), coincides with 
the direction ipCM of the center of mass velocity Vcm, 
i.e., ipi = ipcM at all times. This case can be solved 
exactly for the order parameter 'F, and may arise from 
taking ji = N — 1. The only solution that the system 
presents is that one for the TranS regardless of the initial 
conditions, and can be calculated by solving the equation 
of motion for the center-of-mass. The latter is obtained 
by summing Eqs. fl} and dividing by the number of par- 
ticles N, 



dV, 



rn- 



CM 



dt 

N 

-It 

i=i 



7 v 

V CM 



cos(V>i +<^)i + sin(^ + (Al) 



where the effects of angular noise are captured in the 
expression cos(V>i + <pi)i + sin(-0i + (pi)}, with ipi and (pi 
defined as before. 

We have then, Eq. (|Al[) is akin to the Langevin equa- 
■ where the random force 



tion m 



dt 



£ is defined as the random quantity Yli=i [cos(ipi + 
4>i)i + sm(ipi + </>i)j], and can be interpreted as a di- 
rected random walk with persistence [301 ] . Equation 
(|Aip may be split into two coupled equations for the 
magnitude Vcm = \Vcm\ an d the direction ipcM 
arc: ;u: Ylf V v,i/Y^i v x,i of the center-of-mass velocity. 



After taking an ensemble averages over different realiza- 
tions of the stochastic variables </>; , the equations may be 
explicitly written as 



d(V c 



CM) 
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(Vc 



CM) 



7 

mN 
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COS ( 



(A2a) 



and 



VcM ^r) 



7 
mN 



N 

»=i 



(A2b) 

with a PDF P'(0) 



respectively, where 0i = i/>i — 4>cm 
equivalent for all of the particles. 

In particular, for the global alignment interaction case, 
equations (|A2|) may be solved straightforwardly since 
ipi = ipCM at all times, and 4>i = 4>i (thus P'(4>) — -P»)(0), 
the latter corresponding to the distribution of the noise) . 
Then, the rightmost terms reduce to (cos cf>i) = (cos 4>i) 
and (sin0i) = (sin 4n). 

By definition the probability distribution of the noise, 
P I) (0), given in Eq. ([5]), is the same for all the angles 
(f)i] therefore the quantities (cos and (sin^i) do not 
depend no the sub- index i, and have the values ^ sin 77 

and 0, respectively. The solution to Eq. (|A2al 
given by 



n 

is then 



(V CM (t)) = 
(V C m(0)) e 



— sin 77 

V 



1 



(A3) 



Under the assumption that the system is ergodic, we will 
replace the temporal average in Eq. (|8aj) by an ensemble 



average in the stationary state (Vq M ) which, from ex- 
pression (|A3[) . has the value vo^p- Finally, we get the 
following exact result for the order parameter: 

* = - sin 77. (A4) 
77 

This expression is compared to results obtained from nu- 
merical simulations in the inset of Fig. [Til The agreement 
between theory and numerical data is extremely good, 
which also validates the numerical integration scheme 
used in this work. Notice how this solution only depends 
on the noise intensity 77. 

It is interesting to note that the order parameter in 
the Vicsek et al. model Q, when a global alignment 
interaction is considered in the thermodynamic limit, 
i.e., when N — > 00 for a constant density p = N/L 2 , 
shows exactly the same dependence on 77 as that of Eq. 
(|A4|) . Indeed, in Ref. 0, the instantaneous order parame- 
ter ^(t) = \VcAi(t)\/vQ is also a measure of the coherent 
behavior of the system and by definition given by 
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FIG. 11: (Color online) Normalized magnitude of the mean 
velocity of the group, in the steady state, given by the order 
parameter $ of Eq. (|A5[) , for the Vicsek et al. model, as a 
function of the noise intensity 77 while the limit N — > 00 is 
approached. Numerical (curves with clear symbols) and ana- 
lytical (dashed lines) results, the latter corresponding to the 
expression of Eq. ()A8|) for different values of N, are presented 
assuming a global alignment interaction and a constant den- 
sity p = 10 3 . On the other hand, when a local alignment 
interaction is assumed, numerical results (curves with solid 
symbols) are presented as iV — ^ 00 by making p — * 00 for a 
fixed L = 10. The black solid curve corresponds to the ana- 
lytical approximation of Eq. ()A9|) . In all the simulations for 
the Vicsek et al. model only the small velocity regime [26| was 
considered (vq — 0.03) with r = 1. The inset shows numerical 
results (grey circles) for model JTJ with N = 200, k — 1, and 
t c = 10, when a global alignment interaction (fi/N = N — 1) 
is considered, against the analytical exact result (black solid 
curve) of Eq. ()A4|) . 



where the phases (fi correspond to the directions of the 
individual velocities of N particles contained in a box of 
size L with periodic boundary conditions. The system 
considers local alignment interactions (see below), while 
the velocities Vj of the particles are determined simul- 
taneously at each time step. Then, their positions are 
updated through 



^{t + At) =x i (t) + v i (i)At 



(A6a) 



Here, the velocity was constructed to have a fixed mag- 
nitude vq equal for all of the particles with At =1. In 
this way, the particles advance the same fixed distance in 
the direction of the angle ipi(t + At) at each time step, 
the latter calculated from 



<Pi(t + At) = {<pi(t)) r + <k(t), 



(A6b) 



where {(pi(t)) r denotes the mean direction of the veloci- 
ties of the particles (including particle i) that lie within 
a circle of radius r surrounding a given particle, while <pi 
is a random angle taken from a flat distribution between 
[—77, 77] (here defined in this way to be able to compare 
with our model), with 77 running from to ir 26]. 

It is possible to give an expression for ^> only in terms 
of the random quantities i.e., by assuming a global 
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alignment interaction such that (ifi(t)) r — tp(t), where 
tp(t) is the same for all of the particles and corresponds 
to the mean direction of the whole group. This can be 
achieved by making the range of the alignment interac- 
tion larger or, at least, equal to the size of the system 
(r > L). After taking the sine and cosine of Eq. (lA6bj) . 
and summing over all particles, one gets 



N 



N 



i=l 



sin tfi (t + At) = smip(t) ^ cos <j)j(t) 

i=l 

N 

+ cos ip(t) sin^j (t), 



and 



N 



y^cos ifi(t + At) 



i=l 



COS 1p(t) COS (t) 

i=l 

N 

— s'mip(t) y^ sin 4>i{t). 



By substituting last expressions in Eq. (|A5|) . expanding 
the squares, and simplifying terms one can write 



^{t + At) 
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N 



N 



y^ cos (f>i(t) 



TV 



1/2 



(A7) 



An evaluation of (^(t + At)) in this form is difficult. 
Instead, let us compute the quantity N 2 (^ 2 (t + At)), 

N 2 (V 2 (t + At)) = ^^cos<M^ 



TV 



y^sin cf>i(t) 



Making use of the multinomial theorem [311 ]. after some 
algebra one gets: 

. 2 N-l N 

N 2 (^ 2 ) = N + 2 S -^IJ2J2 1 



= N + N(N - 1) 



i— 1 j>i 

sin 2 r\ 

„2 ' 



and thus 




1 \ sin r] 



N 



(A8) 



In the limit when TV — > oo, and on the basis of some 
observations from our numerical results, one can write 



(A9) 



The last approximation corresponds to the asymptotic 
behavior, in the steady state, for the order parameter of 
the Vicsek et al. model under the assumptions of a 
global alignment interaction in the thermodynamic limit 
(N — > oo). Under the same global alignment assump- 
tion, the expression of Eq. (|A9j) has also been derived by 
approximating the system with a directed random walk 
with persistence [32j . 

Figure [TT1 shows plots of results obtained from numeri- 
cal simulations (curves with clear symbols) for the behav- 
ior of ^>(ri) against the analytical expression of Eq. (|A8|) 
(dashed curves), considering a constant density p = 10 3 , 
and system sizes such that L < r as N tends to infin- 
ity. As apparent from the figure, the order parameter 
approaches the expression of Eq. (|A9|) from above, and 
the absence of the phase transition is evident. For all the 
numerical simulations of the Vicsek et al. model 0, we 
have set At = 1, r = 1, and vo — 0.03, the latter well 
inside the small velocity regime (vq < 0.1) (26j. 

In this regime, particles that were within the same in- 
teraction vicinity at time t will most likely remain within 
the same interaction vicinity at time t + At. Nonetheless, 
in the computation of Eq. (|A9|) we have made no assump- 
tions on the particles speed vq, nor the linear size of the 
system L. Moreover, the same expression has been ob- 
tained for the vectorial network model (VNM) of Aldana 
et al. 33], in Eq. (23b) of Ref.@, when only intrinsic (or 
angular) noise is considered for an infinite network con- 
nectivity. This case corresponds to the global alignment 
interaction limit in both: our model, given in Eq. ©, and 
the Vicsek et al. model. The VNM of Aldana et al., as 
stated there, may be considered as a mean-field approxi- 
mation of the Vicsek et al. model in the limit of extreme 
particle speeds (vq ~ 1000 for simulation purposes) [26| : 
indeed, when vq — > oo. One of the conditions for the 
validity of this mean-field approximation is that spatial 
correlations developed due to the motion of the particles, 
and present in the small velocity regime, are destroyed at 
every time step when the particles are allowed to move 
with extreme speeds. However, as mentioned before, our 
calculation on the Vicsek et al. model does not assume a 
particular regime for Vq. 

The common behavior displayed by ^(77) in these 
three models, is only due to the presence of the global 
alignment interaction in combination with angular noise, 
which in fact are the only features they share. This im- 
plies that the global alignment interaction aids to the 
homogenization of the systems, as it dominates over any 
other characteristic dynamics of a given particular model, 
inducing an effective mean field for all of the interacting 
particles. 

On the other hand, when local alignment interactions 
(r < L) are assumed for the Vicsek et al. model, it is 
well known that the critical point r\ c (that depends on p) 
moves to higher values as the density is increased, finally 
reaching the "infinite temperature" limit, i.e., r\ c — > 7r 
as p — > 00, where the phase transition disappears [3]. 
For this case, and with fixed L = 10, our numerical sim- 
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ulations show that the order parameter approaches the 
expression given in Eq. (|A9|) from below (curves with 
solid symbols in Fig. [TTj) . Notice how the critical point 
rj c of the phase transition is shifted to higher values with 
the increasing p. As density tends to infinity, an thus 
N — ► oo, in the thermodynamic limit, it is plausible 
to expect the fulfillment of the condition for the global 
alignment interaction [(fi(t)) r = ij)(t)] once more, thus, 
leading to the same result as the one presented in Eq. 

(EH). 

This may be explained from the combination of two 



different effects. On the one hand, because of the pres- 
ence of noise in the system, particles move randomly with 
some correlation due to the local alignment interactions. 
As density is increased, this "random" motion makes the 
system spatially more homogeneous, in contrast to the 
low-density case where particles tend to form groups that 
move coherently in random directions [23j | . Then, the 
overlap of the interaction vicinities of the different par- 
ticles increases, and the system becomes more coherent 
developing the effective mean-field interaction previously 
discussed for the global alignment interaction limit. 
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